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Abstract 

Ql^ Room Temperature Ionic Liquids (RTILs) have attracted much attention in the scientific com- 

munity in the past decade due the their novel and highly customizable properties. Nonetheless 
their high viscosities pose serious limitations to the use of RTILs in practical applications. To 
elucidate some of the physical aspects behind transport properties of RTILs, extensive classical 
^1 molecular dynamics (MD) calculations are reported. Bulk viscosities and ionic conductivities of 



butyl-methyl-imidazole based RTILs are presented over a wide range of temperatures. The depen- 
dence of the properties of the liquids on simulation parameters, e.g. system size effects and choice 
of the interaction potential, is analyzed. 
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I. INTRODUCTION 



In recent years, room temperature ionic liquids (RTILs) have see much interest due to 
their promising properties in "gree chemistry" apphcations [jj. Similarly to the well known 
high temperature molten salts, RTILs are liquids composed solely by ions. As opposite 
to molten salts, the presence of large asymmetric organic cations inhibits crystallization 
and allows these salts to be liquid at temperature as low as 100°C. Despite the fact that 
the first report of a RTIL dates back to the beginning of the last century [2], the first 
generation of RTILs, based on chloroaluminate(III) systems, was theoretically predicted |^ 
and experimentally realized {4, 5I just in the 1980's. Nonetheless, practical applications of 
chloroaluminate ionic liquids were strongly limited by their high moisture sensitivity 0, Q] • 
It is just with the advent of water and air stable RTILs [8] that the full potential of these 
new compounds became apparent. Several reviews on this topic appeared in the last years 



m 



clearly highlighting the interest in this field . 

Properties of RTILs are as different from standard molecular solvents as ionic crystals 
differ from molecular ones. Due to the ionic nature of their constituents, RTILs generally 
show negligible vapor pressure and high thermal stability, and they tend to be very good 
solvent for most organic and inorganic species. Moreover, they typically present very high 
electrochemical stability and are intrinsically able to conduct electrical currents. Despite 
all these remarkable properties, none of the existing RTILs possess all of them at the same 
time. Nonetheless, since RTILs are based on organic molecules, they can be easily modified 
by standard chemical reactions. Indeed, both the cation and the anion can be individually 
modified in order to tune the physico-chemical properties of the resulting RTIL. This allows 
a great degree of fiexibility in the design of the most suitable compound, given that the 
number of different possible combinations scales as the product of the number of oppositely 
charged ions. When considering the number of binary and ternary mixtures of all available 
cations and anions, the range of possibilities quickly diverges. In this perspective, the idea 
of task-specific RTILs has emerged and applications of RTILs in the most diverse fields have 
been proposed, and patented. Organic synthesis and catalysis, extraction, and treatment of 
rare-earths elements are some of the most studied fields of application of RTILs. 

Particular attention has been devoted to the possible applications of RTILs in electro- 
chemistry [13\. Specifically, RTILs as electrolytes in dye-sensitized solar cells lUQ, fuel 



cells[18j, and lithium batteries [19jhave been extensively studied both for safety and effi- 
ciency. Indeed, the low volatility of these liquids insures their safety against combustion 
and explosion. On the other hand, some of the existing RTILs present very high electro- 
chemical stability, thus allowing the use of higher voltages in batteries and improving the 
overall efficiency of the devices. A major limitation in actual applications of RTILs as elec- 
trolytes is represented by the high viscosities of all RTILs known to date. Even though the 
conductivity of some of the available RTILs is already sufficiently high for battery appli- 
cations [19], design rules to improve ionic conductivity are currently the subject of intense 
experimental and theoretical study. In this context, much of the scientific effort is directed 
toward understanding of the molecular mechanisms that inffuence the transport properties 
of RTILs. 

Despite the huge increase in the literature on RTILs, experimental results still suffer from 
some limitations, and in particular can be affected by the presence of water and impurities 
in the systems studied [20||. Indeed, water was shown to strongly affect transport properties 
of RTILs, with a decrease of diffusion coefficients by orders of magnitude in water-RTILs 
mixtures being reported [21]. A full quantitative characterization of the amount and kind 
of impurities is often lacking and results can be difficult to compare. 

In order to provide a detailed picture of the molecular mechanisms behind the properties 
of RTILs, several theoretical works have appeared in the literature. The first semi empirical 
and ab-initio calculation on RTILs were performed at the end of the 1980s jsl, 0], with the 
main purpose of investigating the stability and the structure of isolated ions; theoretical 
literature on RTILs started increasing after 2001, when some of the first parametrization of 



22l-l28l| . The first calculations typ- 



classical interaction potentials (force fields) were reported 
ically also reported results on single ions and ion pairs as obtained from density-functional 
theory (DFT). These results were used to fit classical empirical potentials that were exploited 
to compute some of the fundamental properties of the liquids. In order to validate the force 
fields against available experimental results, radial and angular distribution functions, den- 
sities and diffusion coefficients were the first quantities investigated. The agreement between 
simulation results and experimental data were not always satisfactory, with some approxi- 
mation in the functional form of the interaction potential being too crude j29|. In particular, 
it was shown that using an All Atom picture instead of a Unified Atoms Method was more 
effective in reproducing experimental densities ^]. Similarly, including polarization effects 



was shown to significantly improve the agreement with experimental results, in particular 



for transport properties 



30||. Despite some of the limitations of the available force fields, MD 



simulations were effective in predicting and validating some of the new features characteris- 
i;ic of the RTILs' structure, such as the presence of heterogeneities and holes in the liquids 
3l| . Nonetheless, some properties still lack an accurate analysis and characterization. This 
is particularly true for transport properties, such as viscosity and conductivity, that, due to 
the slow glassy-like dynamics of most RTILs, require long simulation times and large system 
sizes. Indeed, while parametrized force fields are generally checked against structural prop- 



32 



erties, few groups have reported careful analysis of viscosities and conductivities 
For these reasons, in the present work an extensive set of classical MD simulations on pro- 
totypical RTILs are reported. After a short methodological section (II), our results on the 
effect of temperature and simulation parameters on the transport properties are presented 
(IIIA). Eventually, results on ionic liquids based on different anions are compared in Section 
IIIB. 



II. SIMULATION DETAILS 



Classical mo' 
throughout [37- 



ecular dynamics (MD), as implemented in the DL_POLY program, was used 



39l |. Three different 3-butyl-l-methylimidazolium (BMIM) based RTILs were 
examined, namely the salts obtained with the PFg , BFj , and bis(trifluorometylsulfonyl)- 
imide (Tf2N~) anions. As for the interaction potential, three different force field were consid- 



ered for the case of BMIM-PFe 



2^, 



40H4a, while simulations on BMIM-BF4 and BMIM-TfsN 



4 CllJ 



were performed using the force fields developed by Canongia-Lopes et al. [4ll-l43|. An inter- 
action cutoff radius of 15 A was used throughout. Initial configurations containing 128, 432 
and 1024 ion pairs were generated starting from a fee cubic lattice with the ions occupying 
random lattice sites. Equilibration in the NPT ensemble was enforced using Berendsen's 
thermostat and barostat j^l, with equilibration runs performed until convergence of the 
statistical average of the density was achieved. Following equilibration, up to four different 
production runs for each system were performed in the NVT ensemble, using the Nose- 
Hoover thermostat [45|, |46| . A simulation timestep of 2 fs was used for temperatures lower 



or equal to 500 K, while 1 fs was chosen for simulations at higher temperatures. 

Diffusion coefficients of the ions were computed in terms of mean square displacements 



(MSD) 

[rHt)) = {\r^{t)-r,{0)\') 



using the Einstein relation |47|, |48 



' t^ooQ dt ^ ^ 



Following standard methodology (e.g. as reported in [3j|), to assess the true diffusive be- 
havior of the ions, MSDs will be displayed in log-log plots, together with their slopes /3 (t) 
as a function of time 

^ dlogt ^ ^ 

At very short times, a value of /3 equal to two is expected, corresponding to a free, ideally 

ballistic, motion of the ions. On the other extreme, at very long times /3 should reach a value 

equal to one, corresponding to real diffusive regime, in which the mean square displacement 

of the ion grows linearly with time. At intermediate times a sub-linear behavior is expected, 

characterized by a logarithmic slope /3{t) lower than one. We would like to stress here that 

this analysis is crucial in determining the effective achievement of diffusive behavior and to 

determine the portion of the MSD plot that should be used to fit diffusion coefficients. This 

is particularly important in the case of slow, viscous liquids such as RTILs. 



The Einstein formalism was also used in the calculation of ionic conductivities 



Mm 



^ = 1™ ^^TITTF— (3) 



t^oo GVkBT dt 
while viscosities {rj) were computed from the auto correlation o: 



of the stress tensor a^^, exploiting the Green-Kubo relation 



the off diagonal elements 



V = ^ I {a^H0)<7^Ht))dt. (4) 



III. RESULTS 
A. BMIM-PF6 

A system composed of 128 BMIM-PF6 ion pairs, described with the interaction potential 
parametrized by Liu et al. (ref. [40], abbreviated in the following as LHW2004), was 
equilibrated at several temperatures, ranging from 300 K to 1000 K. Densities, reported 



Systems 


p(g/cm3) 


D+ (mVs) 


D- (mVs) 


rj (mPa s) 


A(S/cm) 


298.15 K (exp [50]) 


1.37 


1.6 ±0.2 • 10-11 




196 ± 20 


4.04 ± 0.4 • 10-3 


300 K (exp [49]) 


1.37 


8.0 ±3.0- 10-12 


5.9 ±2.4- 10-12 


230 ± 65 


1.66 ±0.3- 10-3 


300 K 


1.36 


4.8 ±1.0- 10-13 


2.6 ±0.5 -10-13 


2400 ± 2000 


1.5 ± 0.7 - 10-^1 


350 K 


1.33 


4.0 ±1.0- 10-12 


2.0 ±0.8 -10-12 


940 ± 400 


1.9 ±0.8- 10-3 


375 K 


1.32 


1.0 ±0.1 • 10-11 


6.2 ±1.0 -10-12 


421 ± 200 


1.4 ±0.5- 10-3 


400 K 


1.29 


3.1 ±0.1 • 10-11 


2.0 ±0.3 -10-11 


111 ±4 


8.7 ±0.6- 10-3 


450 K 


1.26 


1.14 ±0.06 • 10-1° 


7.9 ±0.8 -10-11 


31 ± 1 


1.5 ±0.3- 10-2 


500 K 


1.22 


2.5 ±0.2 • 10-1° 


1.9 ±0.1 - 10-1° 


12 ± 1 


4.8 ± 1.0 - 10-2 


600 K 


1.16 


8.5 ±0.1 • 10-1° 


7.5 ± 0.04 • 10-1° 


7±1 


9.7 ± 1.0 - 10-2 


1000 K 


0.87 


7.6 • 10-*^ 


7.3 - 10-*^ 


5 


2.25-10-1 



Table 1. Densities p, cation diffusion coefficients anion diffusion coefficients viscosities 77, 
and conductivities A of the BMIM-PF6 RTIL as a function of temperature, as computed by classical 
NVT molecular dynamics simulations on a periodic cell of 128 ion pairs using the LHW2004 force 



field. 



in Table [11 are in aefreement with what reported in the literature and show a efood match 

n n 

with the experimental results [49|, [SOfl available at the lowest temperatures. The radial 
distribution functions of the system (reported in Figure [1]) are also in very good agreement 
with previous results [40]. In order to provide a more detailed description of the structure 
of the liquid, when analyzing radial distributions the cation has been subdivided in two 
distinct regions, corresponding to the charged aromatic ring (R) and the long alkylic tail 
(T). Charge-induced correlations in the liquid can be clearly evinced from the reported plots 
and persist, with some broadening, even at the highest temperatures considered. Van der 
Waals interactions, instead, are responsible for the high degree of correlation between the 
alkyl chains of the cations. Due to the local, short-range nature of the Van der Waals 
interaction, the first peak in the tail-tail plot is one of the most sensible to temperature, 
shifting towards larger distances as the simulation temperature increases. 
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Figure 1. Radial distribution functions of BMIM-PFg RTIL, as computed by classical NVT molec- 
ular dynamics simulations on a periodic cell of 128 ion pairs described by the LHW2004 force 
field. The BMIM cation (top panel) is partitioned in two distinct subregions, corresponding to the 
charged aromatic ring (R) and the long alkylic tail (T). The radial distribution functions at 300 
K are reported in black, while for higher temperatures the differentials radial distributions with 
respect to 300 K have been reported (colored lines). Strong charge induced correlations are evident 
from the leftmost series of panels. A reasonably strong correlation is also present between the alkyl 
tails of the cation, as shown in the T-T plot (rightmost series, central panel). All the distribution 
functions show non-negligible correlations for distances as long as half the cell size, even at the 
highest temperatures. 

1. Transport Properties 

Mean square displacements of the different ions are reported in Figure [2] as a function of 
time, for the range of temperatures studied. 

Despite the remaining noise in the logarithmic derivatives at long times, it is seen from 
these plots that a diffusive behavior is reached within the simulation times at each temper- 
ature but the lowest one. Diffusion times, i.e. the times necessary to the ions to reach a 
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Time (ps) Time (ps) 

Figure 2. Mean square displacements (left panel) and corresponding values of /3 (Eq. [2j) for the 
BMIM-PFe RTIL, as computed by classical NVT molecular dynamics simulations on a periodic 
cell of 128 ion pairs described by the LHW2004 force field. Cations are represented by thick lines, 
while thin lines are used for the anions. In order to characterize MSDs at short times, a separate 
set of simulations was used, characterized by a very short (0.5ps) time between sampling different 
configurations. 

diffusive regime with /3 — 1, are summarized in Figure [31 The curves show a divergent behav- 
ior at low temperature and can be fitted with good accuracy by a Vogel-Fulcher-Tammann 
relation / = Aexp{B/{T — Tq)). A value of B 1720/1730 and a divergence temperature 
To ^ 200/215 K can be extracted for the diffusion times of the cation and the anion re- 
spectively. By extrapolating these behavior to the lowest temperature considered (300 K), 
a diffusion time of the order of 10^ ps can be estimated. Thus, a straightforward determi- 
nation of the diffusion coefficients of the ions at this temperature appears beyond the limits 
of standard computational resources. For this reason, the different transport properties of 
the system at 300 K are only estimated by averaging the results over several independent 
calculations. In order to improve the statistics in the results and get an estimate of the 
errors, the same approach was used for simulations below 600 K, that show diffusion times 
of the order of few nanoseconds. 

Diffusion coefficients of both the cations and the anions, reported in Figure [4] and in Table 
[H show Arrhenius behavior at high temperatures, while a deviation from linearity is evident 
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Figure 3. Diffusion times, i.e. times necessary to the ions to reach a diffusive regime with /? 1, as 
extracted from classical NVT molecular dynamics simulations on a periodic cell of 128 BMIM-PFg 
ion pairs, as described by the LHW2004 force field (filled circles). Non-linear curve fits of the form 
/ = Aexp(B /(T — To)) are added to the computed results and represented with solid lines. 

at the lowest temperatures. As expected, diffusion coefficients at 300K appear overestimated, 
due to the lack of a true diffusive regime in the simulations at this temperature. Activation 
energies of ^ = 31.4/28.4 ■ IQ-^kJ/mol for the diffusion of the cation/ anion can be inferred 
from the high temperature slopes of the curves reported in Figure [H The change in slope 
at low temperature is an indication of the glassy behavior of the liquid, consistent with the 
very slow non-ergodic dynamics of the systems in this range of temperatures. As was done 
for the diffusion times, fitting the data with a Vogel-Fulcher-Tammann relation points to a 
divergence temperature Tq ^ 230/235 K. The correct qualitative behavior of the results is 
well reproduced, together with accurate slopes as a functionn of temperature, with anions 
diffusing slower than cations, despite their smaller size. The higher conformational flexibility 
of the cation, due to the presence of alkyl chains, is generally acknowledged as the main 
explanation for this behavior. Nonetheless, computed results are one order of mamitude 

n 

far from the experimental data available at the lowest temperatures [49||. This discrepancy 
may be due to the intrinsic deficiency of the non polarizable force-field to reproduce the 
real dynamics of the systems [sol. On the other hand, the presence of even a small fraction 
of impurities in the experimental samples can be responsible for a significant change in the 
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Figure 4. Diffusion coefficients of cations (black) and anions (red) as extracted from the MSDs 
computed on a periodic cell of 128 BMIM-PFe ion^ pairs, as described by the LHW2004 force 
field (filled circles). Experimental results from Ref. [49] are also reported for comparison (empty 
squares) . 



reported results [2l|. 

It is important to note that original results on the LHW2004 force field |4o| showed an 
impressive agreement with experiments. This apparent contradiction can be explaine d by 
the very short simulation times used to extrapolate the diffusion coefficients in Ref. (40|. 
Especially for results at the lowest temperature, according to the considerations above and 
the results in Figure [3l a sub-linear non-diffusive behavior of the MSD should be inferred 
for all simulation times reported in Ref. 4o|. Longer simulations on the same systems at 



350 K appeared in the literature 



5ll | and showed a remarkable agreement with the present 



results, thus validating our conclusions. 

Viscosities and ionic conductivities of the system are reported in Figure [5] and summarized 
in Table [B As for the case of diffusion coefficients, also viscosities and conductivities show 
significant deviations from experimental results. In both cases, results one order of magni- 
tude lower than the values reported in the literature were obtained, even if their trends with 
respect to temperature are also well reproduced. It should be stressed that the uncertainties 
in the computed values for these quantities are higher than the average errors in computing 
diffusion coefficients: this is due to the fact that viscosity and conductivity are global prop- 
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Figure 5. Bulk viscosities (left) and ionic conductivities (right) of BMIM-PF6, as computed by 
classical NVT molecular dynamics simulations on a periodic cell of 128 ion pairs described by the 
LHW2004 force field. Experimental results from Ref.|49|(exp^. empty squares) and Ref.js^ (exp^, 
empty diamonds) are also reported for comparison. 



erties of the system, while diffusion coefficients calculations benefit from ensemble averages 
over all different particles in the system. 



2. Size effects 

As shown by the radial distribution functions (Figured]), Coulomb-induced ordering of the 
systems is significant for distances comparable and larger than the 128 ion pairs simulation 
cell. Thus, particular care should be taken in evaluating size effects on the structure and on 
the transport properties of the liquid. 

In Figure [6] the radial distribution functions as a function of the number of ion pairs (128, 
432, 1024) are reported at 400 K. Correlations due to the electrostatic interactions in the 
liquid decay completely at half the cell size only for the larger systems, composed by 1024 
ion pairs. Nonetheless, the local structure of the liquid is well reproduced in the smaller 
systems, the differences in the computed radial distribution functions being generally less 
than 2%. 

In agreement with what was found for the distribution functions, results in Table [iT] 
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Figure 6. Radial distribution function of the BMIM-PFe RTIL, as computed by classical NVT 
molecular dynamics simulations at 400 K on periodic cell of 128, 432 and 1024 ion pairs using the 
LHW2004 force field. The BMIM cation is divided in two distinct subregions, corresponding to the 
charged aromatic ring (R) and the long alkylic tail (T) (see Figured]). While correlations involving 
non charged tails disappear for distances as low as 20 angstroms, the radial distribution functions 
of the larger system (black curves, scaled by a factor of twenty) show that charge induced ordering 
is negligible only for distances beyond 35 A. The difference between the radial distribution function 
computed on the larger system and the ones obtained with smaller cells are reported in green (red) 
for the 128 (432) ion pairs systems. Despite the fact that, in the smaller systems, correlation lengths 
are larger than half of the cell sizes, no significant changes in the radial distribution functions are 
found. 

show that system size affects only marginally the computed densities. On the other hand, 
transport properties tend to show more pronounced changes in going from the smallest to 
the largest systems. At the lowest temperature considered (300 K), these differences can be 
attributed to the very slow non ergodic dynamics of the system and to the lack of a true 
diffusive regime. At the highest temperatures, simulations seems to identify the intermediate 
size system (N=432) as the least mobile one. Nonetheless, size effects have the same order 
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Systems 


r(K) 


p(g/cm3) 


D+ (mVs) 


D- (mVs) 


7] (mPas) 


A(S/cm) 


= 128 


300 


1.36 


4.8 ± 1.0 • 10-13 


2.6 ±0.5 - 10-13 


2400 ± 2000 


1.5 ± 0.7 - 10-"^ 


N = 432 


300 


1.36 


1.3-10-13 


5.0 - 10-1-1 


10000 


2.8 - 10-^ 


= 1024 


300 


1.37 


4.4 • 10-13 


4.0 - 10-13 


1700 


1.5 - 10-1 


AT = 128 


400 


1.29 


3.1 ±0.1 • 10-11 


2.0 ±0.3 - 10-11 


111 ±4 


8.7 ±0.6- 10-3 


A^ = 432 


400 


1.295 


2.4-10-11 


1.5 - 10-11 


150 


4.2 - 10-3 


AT = 1024 


400 


1.297 


3.0-10-11 


2.0-10-11 


135 


4.7-10-3 


AT = 128 


500 


1.22 


2.5 ±0.2 - 10-1° 


1.9 ±0.1 - 10-1° 


12 ± 1 


4.8± 1 - 10-2 


A^ = 432 


500 


1.222 


2.5 - 10-1° 


1.9-10-1° 


73 


2.4 - 10-2 


AT = 1024 


500 


1.223 


2.9 - 10-1° 


2.3 - 10-1° 


37 


4.2 - 10-2 



Table 11. Densities p, cation diffusion coefficients anion diffusion coefficients D~ ^ viscosities 
77, and conductivities A of the BMIM-PF6 RTIL as a function of temperature and system size, as 
computed by classical NVT molecular dynamics simulations using the LHW2004 force field. 

of magnitude as the uncertainty in the computed results and values obtained already on the 
smaller systems give reasonable estimates of transport properties. 



3. Effect of the force field 

In order to characterize the effect of the computational parameters on thecalculated 
transport properties of the system, a comparison between different available force fields was 
performed. In addition to the results discussed in the previous sections, simulations with 



the interaction potentials developed by Lopes et al. (Ref. 4l|, |42|, in the following abbrevi- 
ated in CLDP2004), and by Morrow and Maginn (Ref. |23|, in the following abbreviated in 
MM2002), were performed. The differences between these force fields lie mostly in the mag- 
nitude of the electrostatic charges and in the description of the four-body bonded interaction 
(dihedral angles) between the atoms in the ring and the alkyl chains. 

In addition, a few test simulations on the CLDP2004 force field have been performed 
by changing the descriptions of the intramolecular bonds. A set of simulations with the 
harmonic constants of bonds and angles reduced by half is reported, together with results 
obtained by constraining all the bond lengths of the system to their equilibrium values. 
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Figure 7. Radial distribution function of the BMIM-PFg RTIL, as computed by classical NVT 
molecular dynamics simulations at 400 K on periodic cell of 128 ion pairs with different force 
fields. The BMIM cation is partitioned in two distinct subregions, corresponding to the charged 
aromatic ring (R) and the long alkylic tail (T) (see Figure [1]). For each plot, the regions where 
the radial distribution functions present the largest dependence on the interaction potential have 
been highlighted in the insets. While most of the radial distributions show only minor differences 
due to the choice of the force field, the anion-anion distribution function and the intramolecular 
ring-tail coordination show marked changes in the different simulations. Moreover, only one of 
the considered force fields shows the appearance of a shoulder at short distances in the ring-ring 
distribution function. 



For all the reported force fields, computed densities are close to the experimental results. 
On the contrary, radial distribution functions show a marked dependence on the adopted 
force- field (See Figure [7j). The strongest effect of the force field is on the intramolecular 
ring-tail and the intermolecular anion-anion distribution functions. In the more flexible 
molecules, the alkyl chain is able to come in closer contact with the ring, thus creating 
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Systems 


r(K) 


p(g/cm3) 


D+ (mVs) 


D- (m2/s) 


7] (mPas) 


A(S/cm) 


LHW2004 


300 


1.36 


4.8 ± 1.0 • 10-13 


2.6 ±0.5 • 10-13 


2400 ± 2000 


1.5 ±0.7- 10-1 


MM2002 


300 


1.353 


3.2 • 10-12 


2.1 • 10-12 


2200 


4.1 • lO-'' 


CLDP2004 


300 


1.357 


5.0 • 10-13 


2.1 • 10-13 


7600 


2.3 • 10-5 


LHW2004 


350 


1.33 


4.0 ± 1.0 • 10-12 


2.0 ±0.8- 10-12 


940 ± 400 


1.9 ±0.8 • 10-3 


MM2002 


350 


1.305 


1.3-10-11 


9.4 • 10-12 


77 


3.9 • 10-3 


CLDP2004 


350 


1.318 


4.2 • 10-12 


2.3 • 10-12 


903 


7.6 • 10-^ 


LHW2004 


400 


1.29 


3.1 ±0.1 • 10-11 


2.0 ±0.3- 10-11 


111 ±4 


8.7 ±0.6 -10-3 


MM2002 


400 


1.268 


8.0 • 10-11 


5.7 • 10-11 


26 


9.8 • 10-3 


CLDP2004 


400 


1.275 


2.5 • 10-11 


1.7-10-11 


196 


7.4 • 10-3 


LHW2004 


500 


1.22 


2.5 ±0.2 • 10-1° 


1.9 ±0.1 • 10-1° 


12 ± 1 


4.8 ± 1.0 • 10-2 


MM2002 


500 


1.185 


4.7-10-1° 


4.0 • 10-1° 


10 


6.1 • 10-2 


CLDP2004 


500 


1.200 


2.4 • 10-1° 


2.1 • 10-1° 


12 


2.7-10-2 



Table III. Densities cation diffusion coefficients D^^ anion diffusion coefficients D~ ^ viscosities 
?7, and conductivities A of the BMIM-PFg RTIL as a function of temperature and interaction 
potentials, as computed by classical NVT molecular dynamics simulations on a system composed 
by 128 ion pairs. 

some steric hindrance that interferes with the cation-anion interactions. This is reflected in 
the decrease of the second anion-anion peak, together with the slight increase of the mean 
cation-anion separation. 

The effect of constraining all the bonds in the molecule is negligible compared to the 
other parameters in the interaction potentials. On the contrary, relaxing the force constants 
of bond and angles is responsible for noticeable changes in the local structures of the liquid, 
increasing the overall correlation and allowing the charged residues to come in closer contact, 
at the expense of the non-polar portions of the cations. 

Transport properties reflect what was found on the structure of the liquid. In particular, 
the MM2002 potential shows significantly higher mobility with respect to the other force 
fields considered. This is probably related to the reduced flexibility of the ring-tail dihedral 
angle of the cation, that favors an unfolded configuration for the alkyl chain and broadens 
the cation-anion correlation function. Despite the differences between the LHW2004 and 

15 



Systems 


r(K) 


p(g/cm3) 


D+ (mVs) 


D- (m2/s) 


rj (mPas) 


A(S/cm) 


CLDP2004 


300 


1.357 


5.0 • 10-13 


2.1 • 10-13 


7600 


2.3 - 10-^ 


CLDP2004 Rigid 


300 


1.358 


2.8 • 10-13 


2.3 • 10-13 


8600 


1.2 - 10-"^ 


CLDP2004 Flexible 


300 


1.346 


7.0 • 10-13 


6.0 • 10-13 


600 


1.4 - 10-"^ 


CLDP2004 


350 


1.318 


4.2 • 10-12 


2.3 • 10-12 


903 


7.6 - 10-"^ 


CLDP2004 Rigid 


350 


1.318 


3.7-10-12 


1.9-10-12 


474 


7.5 - 10-4 


CLDP2004 Flexible 


350 


1.314 


5.1 • 10-12 


2.0 • 10-12 


209 


1.6 - 10-3 


CLDP2004 


400 


1.275 


2.5 • 10-11 


1.7-10-11 


196 


7.4 - 10-3 


CLDP2004 Rigid 


400 


1.280 


3.3 • 10-11 


2.0 - 10-11 


46 


6.5 - 10-3 


CLDP2004 Flexible 


400 


1.274 


3.5 • 10-11 


2.3 - 10-11 


69 


9.2 - 10-3 


CLDP2004 


500 


1.200 


2.4-10-1° 


2.1 - 10-1° 


12 


2.7-10-2 


CLDP2004 Rigid 


500 


1.203 


2.6 • 10-1° 


2.0-10-1° 


28 


2.7-10-2 


CLDP2004 Flexible 


500 


1.196 


2.5 • 10-1° 


2.0 - 10-1° 


22 


5.0 - 10-2 



Table IV. Densities p, cation diffusion coefficients ^ anion diffusion coefficients D~ ^ viscosities 
77, and conductivities A of the BMIM-PFg RTIL as a function of temperature, as computed by 
classical NVT molecular dynamics simulations on a system composed by 128 ion pairs and described 
with the CLDP2004 interaction potential. Different choices of the strength of bonds and angles 



41 



42|; b) an 



have been adopted, corresponding to: a) the one reported in the original paper 
interaction potential with all the bonds constrained to their equilibrium positions (CLDP2004 
Rigid); c) an interaction potential with all the harmonic constants of bonds and angles reduced by 
half (CLDP2004 Flexible). 



the CLDP2004 force fields, diffusion coefficients for the two cases are in close agreement and 
almost one order of magnitude larger than the experiments. Viscosities and ionic conductiv- 
ities show a less clear trend with respect to the force field, probably due to the uncertainty 
in the computed results due to a lack of a proper ensemble averages. Nonetheless, the 
qualitative picture is consistent with the results of the diffusion coefficients, with MM2002 
showing lower viscosities and higher conductivities at all temperatures with respect to the 
other two force fields. Artificially reducing or increasing the rigidity of the bonds and angles 
of the CLDP2004 force field has a limited impact on transport properties, the overall effect 
being within the numerical accuracy of the reported quantities. Contrary to what was found 



16 



above for the different classes of force- fields, increasing the overall flexibility of the cation 
seems to increase the overall diffusivities in the system. 



B. Effect of the anion 



Three of the most studied anions in RTILs have been considered, namely the PFg , BF4 
and Tf2N~ anions. In order to treat all the systems on similar footings, force fields developed 



by the same group of authors (Lopes et al., Refs. |41l,|42|| and j43|) were chosen for all the ions. 
In examining the radial distribution functions, we can see that the two inorganic anions show 
very similar behavior, with almost identical first peak positions. This reflects the comparable 
size of the two anions and the similar interactions in the two systems between the fluorine 
atoms of the anions and the cation. On the contrary, the larger and more asymmetric Tf2N~ 
anion presents correlations between charged residues that are broader and shifted towards 
larger distances. As expected, the cation-cation distribution functions are the least affected 
by the choice of anion, in particular in the short-ranged tail-tail interactions. 

Densities are in good qualitative agreement with experimental results at room tempera- 
ture, but quantitative results are obtained only for the BMIM-PFg system. The other two 
salts show deviations from experiments much larger than the experimental uncertainty in 
the data. Similarly to what was obtained for the BMIM-PFg RTIL, transport properties 
of the liquids are within one order of magnitude from corresponding experimental data. 
Nonetheless, also in this case, computed results show a general good qualitative agreement 
with experiments at room temperature. The Tf2N salt is the one showing the highest con- 
ductivity and diffusivity and lowest viscosity at room temperature. The PFq anion instead 
is the one giving rise to the least diffusive behavior among the liquids in the whole temper- 
ature range considered. It is interesting to notice that, although computed density changes 
are uniform for the three liquids and linear with temperature, transport properties show 
temperature dependencies strongly related to the choice of the anion. In particular, the 
Tf2N anion is the one showing the smallest variations in the considered temperature inter- 
val. As a consequence, already at 400 K the BF4 liquid shows viscosities and conductivities 
comparable to the organic anion. 
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Systems 


T(K) 


p(g/cm3) 


D+ (mVs) 


D- (m2/s) 


T] (mPas) 


A(S/cm) 


FFq 


300 


1.36 


5.0 • 10-13 


2.1 - 10-13 


7600 


2.3 • 10-5 


BFJ 


300 


1.28 


1.3-10-12 


9.2 - 10-13 


414 


6.1 - 10-^ 


TfsN- 


300 


1.50 


2.5 • 10-12 


1.6 - 10-12 


158 


1.8 • 10-3 


PFg 


350 


1.32 


4.2 • 10-12 


2.3 - 10-12 


903 


7.6 - 10-^ 


BFJ 


350 


1.24 


1.5-10-11 


1.4-10-11 


157 


1.7-10-3 


TfsN- 


350 


1.45 


2.3-10-11 


1.5 - 10-11 


99 


2.2 - 10-3 


PFe 


400 


1.28 


2.5 - 10-11 


1.7-10-11 


196 


7.4 - 10-3 


BFJ 


400 


1.20 


8.5 • 10-11 


6.8 - 10-11 


39 


8.4 - 10-3 


TfaN- 


400 


1.40 


9.7-10-11 


7.3 - 10-11 


40 


1.1 - 10-2 


PFe 


500 


1.20 


2.4 - 10-1° 


2.1 - 10-1° 


12 


2.7-10-2 


BFJ 


500 


1.13 


4.7-10-1° 


4.0 - 10-1° 


25 


7.6 - 10-2 


TfsN- 


500 


1.32 


3.6 - 10-10 


3.2 - 10-1° 


24 


4.1 - 10-2 



Table V. Densities cation diffusion coefficients D^^ anion diffusion coefficients D~ ^ viscosities 77, 
and conductivities A of BMIM-based RTILs as a function of temperature, as computed by classical 
NVT molecular dynamics simulations on a system composed by 128 ion pairs and described with 
the CLDP2004 interaction potential. 

IV. CONCLUSIONS 

In summary, a series of extensive classical MD simulations on BMIM based ionic liquids 
have been reported. A careful analysis of simulation parameters have been performed to 
estimate the accuracy of the computed results. The highly viscous nature of the ionic 
liquids studied poses serious problems in evaluating transport properties and was show 
to require simulation times larger than 100 ns to converge below 400 K. In most cases, 
converged quantities could only be estimated at room temperature, for which the systems 
studied are in fact trapped in glass-like dynamics. Diffusion coefficients, viscosities and 
conductivities all show good qualitative agreement and correct temperature trends when 
compared with experimental results. The correct trend of the transport properties for the 
different ionic liquids is recovered, with BMIM-Tf2N being the less viscous RTIL and BMIM- 
PFq showing the slowest dynamics. Nonetheless, results remain one order of magnitude away 
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from experimental data. This discrepancy could be due to the neglect of polarizability in 
the interaction potentials, or due to even small fraction of impurities in the experimental 
setups. 
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Figure 8. Radial distribution function of the BMIM-PFe RTIL, as computed by classical NVT 
molecular dynamics simulations at 400 K on periodic cell of 128 ion pairs described by the 
CLDP2004 force field. For each plot, the regions where the radial distribution functions present 
the largest dependence on the interaction potential have been highlighted in the insets. Different 
choices of the strength of bonds and angles have been adopted, corresponding to the one reported 
in the original paper 
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42l | (CLDP2004 Original), to an interaction potential with all the bonds 
constrained to their equilibrium positions (CLDP2004 Rigid) and to an interaction potential with 
all the harmonic constants of bonds and angles reduced by half (CLDP Flexible). While the de- 
scription of the intramolecular bonds appear to affect only marginally the structure of the liquid, 
increasing the flexibility of the harmonic three body interactions (angles) has a marked effect on 
the computed distribution functions. In particular, the resulting liquid show a higher correlation 
between the charged residues at the expenses of the short-ranged tail tail distribution. 
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Figure 9. Radial distribution function of BMIM-based RTILs for three different anions as computed 
by classical NVT molecular dynamics simulations at 400 K on periodic cell of 128 ion pairs described 
by the CLDP2004 force field. The three anions considered (A) are PF^ (thick black line), BFJ 
(red line) and Tf2N~(thin green line). The two inorganic anions show very similar distribution 
functions, with almost identical first peak positions. This reflects the similar interactions in the 
two systems between the fluorine atoms of the anions and the cation. On the contrary, the larger 
and more asymmetric Tf2N~ anion presents correlations between charged residues that are broader 
and shifted towards larger distances, while the short-ranged tail-tail interactions in the cation are 
only slightly affected. 
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Figure 10. Density, diffusion coefficients for anions (open symbols-dashed lines) and cations (ffiled 
symbols-continuous lines), viscosity, and conductivity of BMIM-based ionic liquids for three differ- 
ent anions as computed by classical NVT molecular dynamics simulations at 400 K on a system 
composed of 128 ion pairs, as described by the CLDP2004 force field. The three anions considered 



are P 
Ref. 



L_g" (black circles), BF4 (red squares) and Tf2N (green diamonds). Experimental results from 
49I are also reported for comparison. 
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